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ABSTRACT 

Precise subtraction of foreground sources is crucial for detecting and esti- 
mating 21 cm HI signals from the Epoch of Reionization (EoR). We quantify 
how imperfect point source subtraction due to limitations of the measurement 
dataset yields structured residual signal in the dataset. We use the Cramer- Rao 
lower bound, as a metric for quantifying the precision with which a parameter 
may be measured, to estimate the residual signal in a visibility dataset due to 
imperfect point source subtraction. We then propagate these residuals into two 
metrics of interest for 21 cm EoR experiments - the angular power spectrum and 
two-dimensional power spectrum - using a combination of full analytic covariant 
derivation, analytic variant derivation, and covariant Monte Carlo simulations. 
This methodology differs from previous work in two ways: (1) it uses informa- 
tion theory to set the point source position error, rather than assuming a global 
root-mean-square error, and (2) it describes a method for propagating the errors 
analytically, thereby obtaining the full correlation structure of the power spec- 
tra. The methods are applied to two upcoming low-frequency instruments that 
are proposing to perform statistical EoR experiments: the Murchison Widefield 
Array (MWA) and the Precision Array for Probing the Epoch of Reionization 
(PAPER). In addition to the actual antenna configurations, we apply the meth- 
ods to minimally-redundant and maximally-redundant configurations. We find 
that for peeling sources above 1 Jy, the amplitude of the residual signal, and 
its variance, will be smaller than the contribution from thermal noise for the 
observing parameters proposed for upcoming EoR experiments, and that opti- 
mal subtraction of bright point sources will not be a limiting factor for EoR 
parameter estimation. We then use the formalism to provide an ab initio an- 
alytic derivation motivating the "wedge" feature in the two-dimensional power 
spectrum, complementing previous discussion in the literature. 
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Subject headings: cosmology: early universe — methods: analytical — tech- 
niques: interferometric 



1. Introduction 

An understanding of the systematic observational biases in a dataset is crucial for form- 
ing an unbiased estimate of the Epoch of Reionization (EoR) signal, as measured by low- 
frequency interferometric radio telescopes, because the signal is expected to be weak (orders 
of magnitude below contaminating signals and the thermal noise), and the scientific use of 
the result relies on the detailed Fourier structure of the signal. Foreground contaminants to 
the EoR 21 cm signal include bright point sources, confusing (unresolved) point sources, and 
diffuse galactic synchrotron emission. Studies have been undertaken to explore methods for 
removing unresolved sources and diffuse structure, and the effects of these techniques on the 
signal, but it is typically assumed that the bright source subtraction has been performed per- 



fectly (Bowman et al. 2009; Liu et al. 2009 Liu and Tegmark 2011| and references therein). 
These bright sources are often used to calibrate the instrument (e.g., Murchison Widefield 
Array - MWAQ Precision Array for Probing the Epoch of Reionoization - PAPERQ Low 
Frequency Array - LOFAR^J Long Wavelength Array - LWA^ Lonsdale et al.|[2009| |Tingay 
eTaLpOl^l |Parsons et aTpHOt |Stappers et aLpHT] |Ellingson et"al][2009| . Imperfect source 
subtraction may lead to systematic errors in the dataset used to perform EoR measurement. 
It is this effect we study in this paper. 

There have been previous attempts to quantify the uncertainties introduced into an 



image dataset, and consequently a measure of the EoR power spectrum (Datta et al. 2009 



20101), by position errors in the bright point source subtraction. These studies have assumed 



position errors with a global root-mean-square magnitude, and propagated them to the final 
measurement. They have not studied the expected position errors introduced by imprecise 
measurement of the point sources (i.e., uncertainties in parameter estimation due to imperfect 
data). An analysis that derives the expected parameter uncertainties, and then propagates 
them analytically into the visibility and image datasets, will provide the most reliable es- 
timate of the impact of bright source contamination. In this paper we present a method 



1 http: / / www.mwatelescope.org 
2 ht tp : / / eor . berkeley.edu 
3 http: / / www.lofar.org 
4 http: / /lwa.unm.edu 
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for determining the impact of bright source subtraction from first principles. Our method 
provides a straight-forward analysis, based on the Cramer- Rao lower bounds (CRBs), of the 
theoretical level to which the parameters of a source, and therefore our ability to subtract 
it, can be determined, and propagates these uncertainties into the visibility plane and EoR 
power spectra. 

Bright point sources affect interferometric datasets substantially. In addition to sat- 
uration of signal within a synthesized beam-width of the actual source, the incomplete uv 
coverage of distributed synthesis instruments leads to strong sidelobe contamination through- 
out the dataset. For measurement of the weak (~20-50 mK) and structured EoR signal, 
bright point sources (>1 Jy) need to be accurately and precisely removed from the dataset 
before analysis can proceed. Precise estimates of the source parameters are achievable due 
to the strong signal from these sources (high signal-to-noise ratio). Conversely, inaccurate 
subtraction of the brightest sources may lead to large amplitude residuals. The balance 
between these factors will determine the utility of these datasets for EoR science. 

At low frequencies (<200 MHz), the ionosphere produces a measurable perturbative 
effect on the wavefront shape of a propagating signal, even for short baseline arrays specifi- 
cally targetting the EoR (e.g., MWA, PAPER). The ionosphere acts as a phase screen. The 
first order gradient term yields a shift in source position, while higher-order terms cause 
lensing. Spatially compact arrays see very little curvature, making a gradient phase screen 



a reasonable approximation under normal ionospheric conditions. Cohen and Rottgering 



(2009) studied the differential refraction (position change relative to other sources in the 
field) for Very Large Array (VLA) fields at 74 MHz. They found an empirical power law fit 
to the differential refraction as a function of source sky separation of ~0.5 during nighttime 
observations (~0.7 during the day), under normal ionospheric conditions (refraction regime), 
yielding a source shift of ~100 arcseconds at a separation of 25 degrees. In addition to this, 
the whole field suffers an overall shift of ~1— 10 arcminutes. At 150 MHz, one would expect 
the refraction to be half the size, owing to the frequency dependence of phase rotation. 

Wide-field low frequency instruments, such as the MWA, PAPER and LWA, aim to use 
known bright point sources as calibrators to model the instantaneous effect of the ionosphere 
on the wavefront. In addition to these calibration sources, lower flux sources will need to be 
removed from the dataset for EoR signal estimation. In this work we consider an observa- 
tional scheme where the calibration is performed in real-time, necessitating a measurement 
of the instrument and sky response on short timescales, and applying these measurements 
in real-time to the measured data (visibilities). We refer to peeled sources generically as 
"calibrators", although some sources may be peeled as contaminating foregrounds, but not 
used for instrument calibration. 



-4- 



Mitchell et aT| (|2008D list the steps in the MWA Calibration Measurement Loop (CML) 



within the Real-time System (RTS). The CML is designed to track the real-time wave- 
front distortions due to the ionosphere, by performing instrumental calibration on short 
timescales (~8-10 seconds). The calibration is performed on bright calibrators, which are 
then subtracted ("peeled") from the visibility dataset before further data processing. The 
potential for ionospheric changes over the calibration timescale effectively corresponds to a 
re-measurement of the parameters of each point source each 8-10 seconds. Accurate and 
precise subtraction of these calibrators is therefore limited by the amount of information 
available about the parameters of a source (brightness and sky position) in that period, 
yielding a lower limit to the precision with which estimation (and therefore subtraction) can 
occur. For an unbiased estimator, this fundamental limit is set by the Cramer-Rao bound 
(CRB) on parameter estimation (Kay 1993). We will use the CRB to quantify the maximal 
precision with which the parameters (sky position and flux density) of a source may be esti- 
mated, and use these uncertainties to model the residual point source signal in the visibility 
dataset. We then describe a formalism for propagating these errors to power spectra, using 
a fully-covariant analytic derivation. We use this formalism, along with Monte-Carlo simu- 
lations to estimate the residual power and uncertainty in the angular power spectrum and 
two-dimensional power spectrum due to point source subtraction. Note that this approach 
differs from that of Datta et al. (2010) in two ways: (1) it uses information theory to set 
the point source position error, rather than assuming a global error; (2) it propagates the 
full covariance matrix to the power spectra estimates, yielding both the uncertainties on the 
power values, and the level of correlation between fc-modes. We work from first principles 
in our approach. 



Methods 



2.1. Estimation performance: Cramer-Rao lower bound 



To determine the theoretical optimal estimation performance with a given dataset, we 
can calculate the Cramer-Rao lower bound (CRB) on the precision of parameter estimates. 
The CRB calculates the precision with which a minimum- variance unbiased estimator could 
estimate a parameter value, using the information content of the dataset. It is computed as 
the square-root of the corresponding diagonal element of the inverse of the Fisher information 
matrix (FIM). The (ij)th entry of the FIM for a vector of unknown parameters is given 
by: 

"<9 2 logL(x;6> 
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where L denotes the likelihood function describing the likelihood of measuring a dataset, x 
for a given parameter set, 0. For a general covariance matrix, C, N independent samples 
and complex data modelled with signal s, the general expression for the Fisher Information 
Matrix is: 
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for white Gaussian noise with variance, <j (Kay 1993). 



The CRB is a useful metric because it places a fundamental lower limit on the measure- 
ment precision of any parameter. In this work it will be used to gain an understanding of 
the fundamental limits of point source subtraction, and how these impact EoR estimation. 
We have previously used it to quantify the precision with which parameters of a point source 
may be estimated with an interfero metric visibility dataset ( Trott et al.||2011 ). Note that the 
results in this work assume an unbiased and efficient estimator exists, in which case the CRB 
yields the ideal precision. In practise, such an estimator may not exist, and the bound will 
be unattainable. Precision better than the CRB can be achieved with a biased estimator. 



2.2. Calibration field: point source population 

To model a realistic distribution of calibrators, we use the differential source counts of 



Hales et al. (1988), obtained from the 6C survey at 151 MHz: 

^ = aeoo^jy-W- 1 , (4) 

A minimum flux density of 1 Jy will be assumed. We generate a catalogue of flat-spectrum 
calibrators according to equation [4], with sky positions assigned randomly across the field, 
yielding ~400 calibration sources in a 30x30 degree 2 field (~1600 in a 60x60 field), with a 
maximum flux density of 40-60 Jy. 
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2.3. Instrument design 



We focus our attention on the MWA instrument, currently under construction in West- 
ern Australia, with science commissioning expected to commence in the second half of 2012. 
The MWA is composed of 128 electronically-steerable tiles, each of which comprises a ground 



plane with sixteen fixed dipole antennas (Lonsdale et al. 2009 Mitchell et al. 2008 Tingay 



et al. 2012). The maximum baseline will be ~3 kilometres, with an instantaneous observing 



bandwidth of 30.72 MHz over the 80-300 MHz range. The MWA will boast a wide field- 
of-view (~30 degrees at 150 MHz) and high sensitivity, making it an ideal instrument for 
detection and measurement of the 21 cm Epoch of Reionization signal at z ~ 6-11. Table [T] 
lists the fiducial MWA parameters used in this work (<xrms is the theoretical noise level per 
visibility, using the whole instrumental bandwidth). We use the MWA array configuration 



of Beardsley et al. (2012, submitted) for the antenna layout. 



In addition to the proposed tile layout of the MWA, which is composed of a 1.5 km di- 
ameter core of 112 tiles and an outer region of 16 tiles (to the full 3 km maximum baseline), 
we will model a hypothetical instrument with an identical number of antennas and maximum 
baseline, but designed to provide uniform ww-coverage (approximately a Reuleaux triangle). 
We will also model PAPER, another low frequency instrument proposing to perform a sta- 
tistical measurement of the EoR. PAPER currently comprises a 64-dipole instrument in a 
minimally-redundant configuration in the Karoo region of South Africa, with a maximum 
baseline of ~260 m. As well as the current design, Parsons et al. (2011) discussed the ad- 



vantages of building maximally-redundant instruments for EoR estimation, and one of their 
test designs will also be investigated in this work; a square 8x8 grid with antenna separation 



of 26.3 m (yielding the same maximum baseline for both configurations). Figure 2.3 displays 
the antenna configurations for the four arrays considered, and example instantaneous uv 
coverage at zenith. For all experiments we simulate an EoR field at a declination of —30°, 
observing the field over a six hour window per day (hour angle range [-3,3] hours). 



2.4. EoR statistical estimation: angular power spectra 

The primary tools of a statistical measurement of the EoR are power spectra: angu- 
lar power spectrum, two-dimensional [2D] power spectrum, and spherically-averaged power 
spectrum. These metrics quantify the signal power on a given angular or line-of-sight scale 
by forming a coherent + incoherent averaging of visibilities. We will consider the angular 
power spectrum and 2D power spectrum in this work. 



The angular power spectrum is defined by (Morales and Hewitt 



2004 



McQuinn et al. 
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Fig. 1. — (Left) Antenna configurations for the MWA and PAPER considered in this work, 
and (right) instantaneous uv coverage at zenith. 
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2006| |Bowman et"aLl|2009| |Datta et al."]|2010| >: 



^ ^ Ny/v | Vut 

V N 

(uv)El 



(5) 



where iV ul , is the number of visibilities contributing to a given (uv) cell, and the sum is over 
the (m>)-cells contributing to that /-mode (I = 27r|u|). The visibilities, V uv , are coherently 
combined (averaged) within each (uv) cell, while the final estimate associating (uv)-cells in 
/-modes is incoherent (visibilities are combined after squaring). Weighting by the number of 
contributing visibilities forms the maximum-likelihood estimate. 



The 2D power spectrum is given by: 
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where the Fourier modes are in units of inverse distance (Mpc 1 ), and are given by (Peebles 



1993 Morales and Hewitt 2004, and references therein): 

2n\u\ 
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where D M (z), H , f 2 %, z are the transverse comoving distance, Hubble constant, frequency of 
the hyperfine transition (~1420 MHz), and observation redshift, respectively. E(z) is (Hogg 



2000): 



E(z) = y/Sl M (l + z) 2 + + z) 2 + ^a, 
and we assume a concordance cosmology (H =70A kms~ 1 Mpc~ 1 , Qm=0.27, flk=0, fl/ 



(9) 
=0.73, 



Komatsu et al. 2011), and an observation redshift of z=8. 



In practise, the power spectra are formed from an image cube of real quantities. In this 
work, we neglect the effects of gridding and image formation, and form the power spectra 
directly from the measured visibilities. This allows a clean propagation of errors from the 
measured data (i.e., the visibilities) to the final statistical measurement (i.e., the power 
spectrum). Because this involves (at most) one Fourier transform (FT, v — > i] oc fen), we 
need to be careful to simulate the usual three-dimensional FT ([9 X , 9 y , v] — > [u,v,rj\) from 
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image space. A three-dimensional FT from real-valued to complex-valued quantities (as 
occurs transforming from the image plane) transforms from real to complex values. We 
typically choose the frequency direction as the final FT. In general, a real-to-complex FT of 
dimensions NxNyNz transform to a complex cube of dimensions NxNyxNz/2. The final 
dimension contains Nz/2 redundant channels, in the same way as a ID real-to-complex FT. 
The choice of the redundant axis is arbitrary, but for convenience we choose the frequency 
direction. Under this scheme, both visibilities and their conjugates (Hermitian reflection in 
iw-plane) must be included in the power spectrum estimation, and only the non-redundant 
half of the rj array is used. 

Forming the power yields a non-zero expected noise power (oc a v ). As discussed by 



McQuinn et al. (2006) and Bowman et al. (2009), the expected thermal power for proposed 



experiments will exceed the expected 21 cm EoR signal power. However, because the noise is 
stochastic, the expected power can be estimated and subtracted, or the experiment divided 
into two and cross-correlated, to eliminate the signal. These techniques cannot remove all 
of the noise power, because they are limited by sample variance, and power will remain at 
a lower level that corresponds to the uncertainty in the magnitude of the noise component. 
The cross-correlation power spectrum is given by ( McQuinn et al.||2006 Bowman et al.|[2009 



Datta et al. 2010): 



a 



(uv)el 



(10) 



^ y 2_/V u , u 
(uv)el 

where V%, Vi are the coherently-summed visibilities from each dataset. The noise variance 
(sample variance) remains after the expected power is removed. This process takes a coherent 
addition and replaces it with a partly incoherent addition, yielding a factor of two increase 
in the noise uncertainty. 

Similar techniques could also be used to reduce the impact of the residual point source 
signal, assuming the residuals are uncorrelated. If there is a substantial degree of correlated 
noise between different angular modes, this technique will have limited use because it cannot 
remove correlated power, and structure will persist. Thus, it is the noise uncertainty that 
is important for comparison. We demonstrate later that imprecise subtraction does yield 
correlated signal in both angular and spectral space, and this technique is therefore of limited 



use. 
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3. Analysis 

To determine the magnitude of residual signal from point source subtraction in the 
power spectra, we use both an analytic propagation of errors, and Monte-Carlo simulations. 
Ideally, would like to determine the full covariant uncertainties for all results, but the com- 
putational demands of forming large, complex-valued covariance matrices necessitate the use 
of simulations for the two-dimensional power spectrum. 

In this section we begin by describing the construction of the power spectra calculations 



(3.1). We then derive the source position uncertainties based on the CRB analysis (3.2) 



and the covariant propagation of these uncertainties to visibilities (section 3.3). These are 



then used as starting points for both the analytic error propagation (3.4, used for the one- 



dimensional power spectrum), and for the simulations (3.5, used for the two-dimensional 
power spectrum). 



3.1. Construction of estimates 

In a real EoR experiment, the iti>-gridding scale (cell size) is the region over which the 
visibilities can be considered coherent, and is determined by the instrument field-of-view and 
desired sampling of the synthesized beam (yielding a pixel scale in image space). For the 
MWA, and sampling the beam with five pixels, this corresponds to ~1000xl000 uv cells. 
The power in each of these cells is then combined incoherently to form the power spectrum 
estimate for a given angular scale. It is computationally infeasible to perform a covariant 
analysis with a complex matrix with ~10 12 entries. For the purposes of this work, where 
we seek (1) to determine the magnitude of the residual point source signal relative to the 
thermal noise component, and (2) to determine the level of correlation between angular 
modes via a fully covariant analysis, we use a very coarse sampling of the uv plane to form 
the intermediate, coherent uv cells (a 40x40 grid of cells, yielding a 1600x1600 complex 
covariance matrix) for our covariant analytic computation. This yields an estimate of the 
angular power that is more coherent than is possible with a real experiment, effectively 
coherently combining visibilities that are not coherent in practise. In addition, we coarsely 
bin the angular scales (/-modes) for computational ease. 

Appendix A contains a derivation of the expected power and power variance for a given 
coherent /incoherent combination of visibilities. This is useful both for an understanding of 
how to design a real EoR experiment, and also in the context of the computational sim- 
plifications we have employed in the analytic covariant calculation. A further consequence 
of this simplification is that the estimated magnitude of the uncertainty (square-root of 
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the sample variance) of the thermal noise power, and of the residual point source subtrac- 
tion signal power will be comparable to the expected power; this directly follows from the 
unphysically-coherent calculation we are performing. 

In addition to the covariant angular power spectrum calculations, we perform an analytic 
variant ID calculation with (uv) sampling appropriate for a real experiment. This yields 
realistic results because we expect the angular modes to be uncorrelated due to Earth rotation 
synthesis over the contributing baselines (structure will be suppressed by the rotation of the 
array baselines with respect to the observation field). This analysis is not possible for the 
2D angular power spectrum, where frequency-based correlations are known to be important 
for reproducing realistic results, and indeed are the basis for the "wedge" feature observed in 
simulations in previously published results. For the 2D power spectrum, we use Monte-Carlo 
simulation. 



The signal in each visibility is the linear combination of signals from each calibrator, 
and is given by: 



where N c is the total number of calibrators, described by source strength, Bi, located at sky 
position, (/j, mi), for a baseline, n, and frequency channel, / (we have assumed flat-spectrum 
calibrators throughout). The signal is embedded within white Gaussian thermal noise, with 
diagonal covariance matrix, C = a 2 I, where a is set by the correlator integration time (dt) 
and channel width (dv), and is given in temperature units by (single polarization): 



3.2. Calculation of uncertainties in calibrator parameters 





a = 



e- 1 sys 



[K], 



(12) 



\ 2 \/dvdt 



where A e is the effective antenna area at wavelength, A, and T sys is the system temperature. 
Using the full dataset, the minimum uncertainty in the parameter estimates for calibrator, 
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and iV is the number of baselines, F is the number of frequency channels. Figure [2] plots 
the optimal point source position precision (Al) as a function of calibrator signal strength 
(Jy) for the four instruments considered. The linear dependence of precision with signal 
strength is observed. The inner core + outer ring structure of the MWA degrades its per- 
formance slightly compared with the uniform array. The short baselines and fewer antennas 
of the PAPER arrays degrade their ability to localise sources well compared with the MWA, 
although their wider instantaneous bandwidth balances this effect somewhat. The vertical 
line at 1 Jy displays the peeling cutoff considered in this work. 

For ionospheric calibration, each calibrator source provides information about the am- 
plitude and direction of the ionospheric phase shift, by comparison with a sky model of the 
true calibrator location. Both the overall and differential refraction are important for a full 
ionospheric model, but the overall shift can be estimated jointly from the calibrators, while 
the local refraction is independent (in the simplest case where local correlations are ignored). 
For faint calibrator sources, the precision theoretically achievable with an 8 second dataset 
may be comparable to the differential refraction due to the ionosphere. In this work, we 
assume no a priori information about source location from previous calibration solutions, 
and use the information contained within the current visibilities alone to estimate source 
precision. This a reasonable assumption for the MWA, but not necessarily for faint sources 
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Fig. 2. — Optimal calibrator position precision (arcseconds) versus source flux density (Jy) 
for the instruments considered in this work, and using the full instantaneous bandwidth of the 
instrument (Table [T) BW = 30.72 MHz (MWA), 70.0 MHz (PAPER); At=8 s, T sys =440 K). 
The vertical line denotes the minimum flux density for peeling (1 Jy). 



-14- 



with PAPER. For this reason, we limit the minimum peeling flux density to 1 Jy in this 



work. In Datta et al. (2010), because all sources are assumed to have the same RMS posi- 
tion precision, the brightest source in the field contributes the greatest residual error. They 
concluded that 0.1 arcsec precision was required. Figure [2] shows that this precision can, in 
principle, be achieved by the MWA with an 8 s snapshot for sources brighter than 30 Jy. 
Sources fainter than this cannot be subtracted with this precision, however the magnitude 
of the error in this case is lower due to the source's lower flux. In the next section we 
demonstrate how sources in the field contribute to the overall residual error. 

For the experiments considered in this work, we use the full bandwidth of the instru- 
ment (Table [TJ to determine the precision of source parameter estimates, regardless of the 
bandwidth used for the individual power spectra. This provides the highest precision on 
parameter estimates available to the instrument. 



3.3. Propagation of errors into interferometric visibilities 

Figure [3] displays schematically the steps involved in forming the power spectrum, and 
in propagating the errors from residual point source signal to the power spectrum. 

The uncertainties in signal parameters for each calibrator are propagated into uncer- 
tainty in the measured visibilities, in addition to the statistical thermal noise. Because the 
measurement errors for each calibrator are independent, the overall visibility-based noise 
covariance matrix will be the sum of the matrices for each calibrator. The error propagation 
follows a fully- covariant treatment, where the residual signal noise covariance for calibrator, 
i, is given by: 

C Vi = JC 9l J\ (20) 

where J is the 3 x N vis complex Jacobian of partial derivatives of the visibility functions with 
respect to the parameters (0j = (!„mj, -Bj)), Cg i is the 3x3 covariance matrix of parameter 
uncertainties, and the dagger denotes the complex-conjugate transpose operation. 

The treatment provides a full covariant analysis, including any correlations between 
visibilities due to residual noise power. For a given visibility and a single calibrator, the 
variance is given by (using equation 20): 
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Fig. 3. — Diagram showing the steps in estimating the power spectrum (upper boxes), and 
errors on the power spectrum (lower boxes). 
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where /i, /2, fz are functions of the baseline lengths (array geometry). The residual error in 
each visibility due to each calibrator is independent of the calibrator signal strength. There- 
fore, the impact on the visibilities of subtracting each calibrator has the same magnitude 
for all calibrators [but different distributions across visibilities; the covariances also depend 
on the source positions, (li,mi)]. This result deviates from the assumptions of previous 
work, which assumed that the source position error was independent of signal strength, and 
therefore that residual error was greatest for the strongest calibrators. 

Because the residual signal is due to a real source (signal-like), and not a random process 
(noise-like), the errors will be correlated between baselines and between frequency channels, 
for a given estimation of the source parameters. These covariances yield a correlated noise 
covariance matrix for the visibilities, deviating from the uncorrelated white noise of the 
stochastic thermal component (diagonal covariance matrix). Such covariances propagate 
from the visibilities to the angular power spectra, potentially yielding correlations between 
angular modes. 



3.4. Propagation of errors into angular power spectra 

As described schematically in Figure [3] and section 2.4[ the high-level steps to calculating 
the angular power spectrum from visibility data are (McQuinn et al. 2006, and references 
therein): 



• Visibilities are averaged into pre-defined wf-cells in each frequency channel, yielding a 
coherent combination of information; 

• For multiple frequency channels (two-dimensional power spectrum), averaged visibili- 
ties are Fourier transformed in frequency to form the line-of-sight power {y— r\ Fourier 
pair); 

• Averaged complex visibilities are squared to form real-valued power estimates in each 
uvrj-cell; 

• The maximum likelihood estimate of the power spectrum is computed by forming a 
visibility- weighted power estimate in each angular and line-of-sight mode (k±, k»), from 
that mode's contributing uvrj cells. 



The final step can be varied in two or three dimensions to yield one- and two-dimensional 
power spectra estimates. 
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The corresponding error propagation from visibilities to angular power spectra follows 
a fully-covariant error propagation, given broadly by: 

• Propagate errors from visibilities to uv cells according to C uv = JCv t J\ where the 
Jacobian contains normalizations for the number of visibilities contributing to each 
cell, and zeroes otherwise; 

• Squaring averaged visibilities to form power estimates yields a quadratic form (i.e., a 
metric quadratic in the data). The Fourier transform from observing frequency (u) to 
line-of-sight frequency (if) can be incorporated into the quadratic form via a weighting 
matrix. The errors for one- and two-dimensional power spectra propagate according 



to flSearle||1971[ ); 

If x ~ CN{0, C) then Ax ~ N[tr{AC), 2tr(AC t AC)] (24) 
where A is a real weighting matrix, and is given by, 

- ID: A = I, 

— 2D: A a p = cos (27r(v a — vp)r)/N), where N is the number of frequency channels. 

• Propagate to the incoherent addition of power samples (the maximum likelihood es- 
timate of the power spectrum) via the error propagation, Cp = JCmJ\ where the 
Jacobian entries contain the number of visibilities contributing to each (k±, k») cell. 

After cross-correlation to remove persistent noise, the covariances calculated according to 



equation 24 form the noise variance that would contaminate the power spectrum (see dis- 



cussion in Section 2.4) 



3.5. Monte-Carlo simulations 

As discussed in section [3j the expected power and power uncertainty for the 2D power 
spectrum is obtained via Monte-Carlo simulation. To produce realistic results that are con- 
sistent with our theoretical framework, we generate our simulated visibilities from the cor- 
related CRB covariance matrix of source parameter uncertainties (including source position 



covariances, equations 13 19). That is, for each antenna configuration, we simulate visibil- 
ities with thermal noise and residual point source signal (obtained by generating Gaussian 
random variables, with a correlation structure given by the CRB results) at the appropriate 
uv co-ordinates for the given experiment. 
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We follow the standard methodology to form power spectra from visibility data (without 
transforming to the image plane), as described in section 2.4 We simulate 100 noise real- 
izations for each experiment, with parameters given by Table [TJ and with (uv) and (k±,k\\) 
sampling appropriate for a real experiment (e.g., 1000x1000 uv cells, and ten k bins per 
decade, for MWA). In this way, we achieve realistic power spectra using the theoretical 
source position errors derived from the bounds. 



4. Results 

We have presented the formalism for deriving the expected power and power covariance 
in angular power spectra. Computing the covariances between angular or line-of-sight modes 
provides an understanding of the level of correlation in the data, introduced by the peeling 
process. In addition to desiring the amplitude of the peeling residuals, compared with ther- 
mal noise, we also require an understanding of the correlations between modes, to perform 
an optimal and unbiased estimate of EoR parameters from power spectra data. As described 
above, performing a fully-covariant analysis is limited by the computational requirements 
of calculating and storing very large matrices. Thus, we also perform variant analyses and 
Monte-Carlo simulations. 

We now present the results of an analytic covariant angular power spectrum, which 
is binned coarsely for computational feasibility and provides information about the level of 
correlation between I modes. In addition to the covariant calculation, we then present an 
analytic variant angular power spectrum, which is binned finely, and reflects the conditions 



available to real experiments (Section 4.1 ). We then present the results of a simulated covari- 
ant 2D power spectrum (where the covariances are carried naturally through the simulation, 
but without the number of simulations required to form an accurate estimate of the mode 
covariances), which is binned finely, and also reflects the conditions of a real experiment 
(Section gg. 



4.1. Angular power spectrum 



Figures [4] and [5] display the fully-covariant, and variant, analytic calculation results 
for the angular power spectrum, for each of the four experiments considered, and with 
experimental parameters defined in Table [TJ 

Figure [4 shows the uncertainty on the residual signal power {^JC~c l ) 1 given by the 
square-root of the covariance term in equation 24 It plots I — I and is symmetric about the 
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Fig. 4. — Square- root of covariance for angular power spectrum (mK 2 ). 
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Fig. 5. — Expected power and uncertainty for thermal noise and residual signal (mK 2 ), using 
a variant calculation with realistic (uv) cells. The binning is logarithmic with ten bins per 
decade in I. The expected contribution from the 21 cm signal is also shown for comparison. 
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diagonal, which represents the variance in each /-mode (note that no factor of 1(1 + 1)/2tv has 
been included here). Covariances are ~2 orders of magnitude below the variances, indicating 
that angular modes are not substantially correlated. This result is unsurprising given that 
we expect that Earth rotation synthesis of the field across six hours of sky will average out 
any correlations between baselines (note that rotation is not expected to remove frequency 
channel-based covariances, which rotate with the other visibilities on the same baseline). This 
result suggests that correlations between angular modes are minimal, and EoR parameters 
derived from measured power spectra will not be strongly biased by assuming uncorrelated 
data. 



5j displays the expected power (1(1 + l)Ci/2n) and power uncertainty 
noise and residual signal, as well as the expected 21 cm signal (assuming a fully- 



Figure 
the therma 

neutral IGM at z=8, where we have used a linear power spectrum and assumed isotropy). 
In all cases the thermal noise power and uncertainty exceed the power from the point source 
subtraction. One can see the differences in the gradients of these noise components, with 
the higher angular modes suffering more severely from point source subtraction, consistent 



with the additional noise in the long baseline visibilities (equation 23). These figures also 
highlight the differences between the array configurations; the minimally-redundant arrays 
(Figure [5^,c) exhibit uniform behaviour across a wide range of angular modes, while the 
non-uniform arrays (Figure [5|3,d) sample particular modes well, but have a reduced range 
of sampled /-modes. In particular, the maximally-redundant PAPER grid array displays 
slightly improved performance at the lowest sampled modes (effectively the grid spacing), 
while performing poorly at the higher modes, and sampling a reduced range of modes. The 
improvement at small I is limited by the logarithmic binning of the /-modes, yielding fewer 
visibilities contributing than would be the case for a linear binning. 

The MWA configurations show similar behaviour: the uniform array yields even per- 
formance across the angular scales, while the centrally-concentrated array samples smaller 
/ modes, yields improved performance at small /, and poorer performance on the longer 
baselines (incorporating the few antennas on the outer ring). 



4.2. 2D power spectrum 

Figures [6^ a-d) and [7] show the uncertainty in power for thermal noise and residual point 
source signal (mK 2 Mpc 3 ), respectively, for the four experiments considered, and obtained 
via Monte-Carlo simulation. Results from the same instrument (MWA or PAPER) are 
binned identically, have the same axis ranges, and the same minimum power level, but color 
scales differ. In addition to the thermal uncertainty, figure [6te) displays the expected signal 
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from the 21 cm emission, assuming a fully-neutral IGM at z—8, where we have used a 



linear power spectrum, and the transfer functions of Eisenstein and Hu (1999), and assumed 



isotropy (Furlanetto et al. 2006 McQuinn et al. 2006). 



The thermal noise (Figure [6]) is uncorrelated between frequency channels, yielding the 
same expected power and uncertainty along the line-of-sight (k\\), where the number of 
contributing visibilities is equal. On angular modes, the uncertainty is set by the number of 
visibilities within a (uv) cell, and the binning used in the plotting (i.e., the coherent versus 
incoherent balance, and the size of the bins). The same behaviour as seen for the angular 
power spectrum is observed: minimally-redundant arrays (Figure |6^,c) yield flat power over a 
wide range of angular modes, while the centrally-concentrated/maximally-reduandant arrays 
(Figure |6}3,d) sample some modes particularly well, but at the expense of overall performance 
and range of sampled modes. 

The residual signal uncertainty (Figure [7]) displays qualitatively different behaviour to 
the thermal noise, owing to the correlated frequency channels and the dependence of noise 
The wedge-shaped feature extending from small k±-smdl\ kn to large k±- 



variance on \u\ 



large kn was first reported by Datta et al. (2010), after initial indications from 



(2009), and has been further discussed by Vedantham et al. (2012); Morales et al. (2012), 



Bowman et al. 



and is the power footprint of the instrument PSF. In Section [5] and Appendix [B] we explore 
the origin of this feature in more detail, and provide a pedagogical derivation of its shape 
and location. 

The residual signal power shows a clear ramping of signal towards smaller angular scales 
(large k±), again resulting from the \u\ 2 dependence of visibility-based noise variance. 



Figure [8] displays the 21 cm signal-to-noise ratio (expected signal divided by noise un- 
certainty) for the actual MWA configuration, and for the thermal noise and residual signal 



contribution. Consistent with previous estimates (Datta et al. 2010 Beardsley et al. 2012, 



submitted), we do not expect to detect the 21 cm EoR signal for 300 hours of observing and 



the parameters used here (SNR<1), when considering thermal noise. There is, however, a 
region of parameter space for which the signal-to-noise ratio exceeds 0.1, suggesting a pos- 
sible detection in this region for ~1000 hours (assuming coherent addition of visibilities). 
For the case of the residual point source signal, there is a large portion of parameter space, 
corresponding to high expected 21 cm signal and/or low uncertainty, where the SNR exceeds 
unity. These results are consistent with discussion isolating the low contamination region of 
k± — k\\ space for EoR detection and estimation. 
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(c) Actual 64-dipole minimum-redundancy (PAPER). (d) Potential 64-dipole maximum-redundancy (PAPER). 




(e) Expected 21 cm HI signal. 



Fig. 6. — (a-d) Thermal noise uncertainty (mK 2 Mpc 3 ). Angular modes are binned logarith- 
mically, with 10 bins per decade, while line-of-sight modes are linear and set by the Fourier 
modes in the Fourier transform. Results from the same instrument (MWA or PAPER) are 
binned identically, have the same axis ranges, and the same minimum power level, (e) Ex- 
pected 21 cm signal assuming a neutral medium at z=8 (see text for assumptions used in 
model). Parameters from Table fl] are used for these results. Note the different color scales. 
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(c) Actual 64-dipole minimum-redundancy (PAPER). (d) Potential 64-dipole maximum-redundancy (PAPER). 



Fig. 7. — Residual signal uncertainty (mK 2 Mpc 3 ). Angular modes are binned logarithmi- 
cally, with 10 bins per decade, while line-of-sight modes are linear and set by the Fourier 
modes in the Fourier transform. Results from the same instrument (MWA or PAPER) are 
binned identically, have the same axis ranges, and the same minimum power level. Parame- 
ters from Table Q] are used for these results. Note the different color scales. 
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Parameter 


Value (MWA) 


Value (PAPER) 


N ant 


128 


64 




150 MHz 


150 MHz 


Bandwidth 


30.72 MHz 


70.0 MHz 


Av w 


125 kHz 


125 kHz 




8 MHz 


8 MHz 


Field-of-view 


30° 


60° 


Calibrators (>1 Jy) 


392 


1579 


T 

sys 


440K 


440K 


At 


8s 


8s 


°RMS 


28.1 mK 


18.6 mK 


^tot 


300h 


720h 


^max 


3000m 


260m 



Table 1: Parameters used for the primary instrument design of the MWA and PAPER 64- 
dipole instruments. For all experiments we use the full instantaneous bandwidth of the 
instrument to estimate source parameter precision. 



Iog, SNR 



log, SNR 





(a) Thermal noise uncertainty, MWA. 



Io9,o (Mpc 

(b) Residual signal uncertainty, MWA. 



Fig. 8. — Signal-to-noise ratio of EoR statistical estimation assuming a fully-neutral IGM 
and isotropy, for (a) expected 21 cm signal to MWA thermal noise uncertainty, and (b) 
expected 21 cm signal to MWA residual signal uncertainty. 
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5. Explaining the "wedge" feature 

In Appendix [B] we provide a full pedagogical motivation for the shape of the "wedge" 
feature in the 2D power spectrum. Here we present the main results, and leave the interested 
reader to follow the full derivation in the Appendix. 



The wedge feature is the footprint of the instrument PSF, transformed into power (Datta 



et al. 2010; Bowman et al. 2009 Vedantham et al. 2012 Morales et al. 2012). More specifi 



cally, it is a combination of two effects, which are then propagated through a Fourier trans- 
form, and squared to form power. The two effects are: (1) the integrated effect on frequency- 
based correlations of the subtraction of multiple point sources distributed randomly across 
the sky; and (2) the spread of frequency channels across (uv) cells for visibilities observed 
with the same baseline (mode- mixing). 

Residual signal from subtraction of a single point source from the visibilities yields 
frequency channels that are fully correlated in amplitude, but with a phase shift that is 
dependent on the source location and the baseline vector. Integrating the effect of multi- 
ple sources yields a sinc-like function in the frequency correlations, effectively reducing the 
frequency-based correlation length. In addition, the spread over coherent (uv) cells of the vis- 
ibilities forming a single baseline (mode-mixing) effectively truncates the sinc-like frequency 
correlations at the first sine null. This effect is most pronounced on the longest baselines, 
and non-existent on the shortest (typically), where the visibilities from a given baseline all 
contribute to the same (uv) cells. The frequency-based covariance matrix is therefore rep- 
resented by a real-valued sine function, multiplied by a rectangular window function, with 
width of the first sine null. The Fourier transform and squaring of this functional form yields 
a wedge shape given approximately by: 

H/(A; ± ,A ;|! )cxsinc 2 f— ^-^-V (25) 

where, 

nc(l + z) 2 D M (z) 

" = 2H f 21 E(z) ' (26) 

are array-independent terms defining the cosmology and observation redshift (with distance 
units), r max is the instrument field-of-view (radial), BW is the experiment bandwidth, and 
v is the centre frequency of the observation. Here W represents the uncertainty in noise 
power (the expected power follows the same functional form) due to imprecise subtraction 
of multiple point sources, spread randomly across the instrument field-of-view. 



We have therefore demonstrated that the wedge shape in the line-of-sight direction can 
be explained by a combination of mode-mixing (truncating the correlation length in fre- 
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quency) and a real-valued sine-like correlation function (derived from addition of random 
phases from individually-correlated point sources at random sky positions). In the angular 
modes, longer baselines are more severely affected than shorter baselines, with an approx- 
imate \u\ 2 dependence of visibility variance (and consequently, expected power and power 
uncertainty). 



6. Discussion and conclusions 

Our results indicate that (1) optimal point source subtraction to 1 Jy will not be a 
limiting factor in statistical EoR estimation for proposed experiments with MWA and PA- 



PER, contrary to the conclusions of previous analyses (Datta et al. 2010), which relied on 
assumed position errors that were independent of source strength; (2) frequency-channel 
covariance combined with spread of visibilities across coherent (uv) cells yields the "wedge" 
feature observed in point source subtracted EoR simulations; (3) angular mode correlations 
are expected to be small (1-2 orders of magnitude below variance level). We have presented 
an end-to-end derivation of the approximate shape of the wedge feature, using the first prin- 
ciples framework we have developed. Here we describe some interesting features and caveats 
of these results, based on the assumptions within our framework. 

The magnitude of the residual point source signal is dependent on the number of calibra- 
tors being peeled, the thermal noise level using the full array and bandwidth, and the angular 
mode in question (with a ~|it| 2 -dependence on variance in the visibility). The uncertainty 
in the power spectrum at angular mode I therefore scales with: 

0PS(O a K M A + A ' 27 

Ai/tot iVvisAtAz/tot 
where cr t hcrm is the thermal noise level across the entire bandwidth, Az/ tot , for power spectrum 
channels of width At/. The magnitude of the expected power, and power uncertainty, there- 
fore scale linearly with the number of peeled sources, and inversely with the total system 
bandwidth and experiment time. The PAPER array suffers from 1/4 the visibilities of the 
MWA and a larger field-of-view, but this is balanced somewhat by the additional bandwidth 
and the smaller I modes available to the (more compact) instrument. 

We have assumed a sequential peeling scheme, whereby sources are independently and 
sequentially estimated, and peeled from the dataset, without downstream impact on the 
estimation of the parameters of other sources. In reality, an experiment can choose to either 
(1) estimate the parameters of all sources simultaneously, whereby the approximation of 
pure radiometric noise in the dataset (ignoring calibration errors) is appropriate, or (2) 
estimate and peel sources sequentially, whereby the covariant noise matrix at each iteration 
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is a combination of uncorrelated radiometric noise and correlated residual signal noise. The 
choice between these estimation techniques depends on the number of sources being peeled 
relative to the available data (number of unknowns versus constraints; the CRB is singular 
for underconstrained problems), and we leave this discussion for future work. If we ignore 
correlations between visibilities, a simple back-of-the-envelope calculation demonstrates that, 
for the parameters used in this work, our assumption is reasonable (where we have effectively 
expanded the noise covariance matrix as a Taylor series, and kept the first linear term). For 
subtraction of a larger number of sources, or for fewer visibilities, this approximation breaks 
down, and the compound effect of imprecisely estimating source parameters within ever- 
noisier data will increase the residual signal power towards the amplitude of the thermal 
noise power. 

A related assumption in our model is the ignorance of previous calibration solutions 
for estimation of the parameters of a source. As demonstrated in Figure [2j optimal source 
position estimates with the full-bandwidth PAPER array (and an 8 second integration) 
have precision comparable to the potential ionospheric differential refraction magnitude, for 
sources close to the 1 Jy peeling floor. In these cases, the utility of the current dataset may 
be augmented by use of prior information about the local behaviour of the ionosphere. For 
larger arrays, this is also an issue if one attempts to peel below ~1 Jy. The optimal method 
for estimating source location (balance of prior and current information, and how these are 
combined) will depend on the peeling floor, the efficiency of the estimation method, and the 
specific array design. 

Discussion of realized estimation performance highlights a further basic assumption of 
our framework; we have proceeded from a fundamental reliance on the existence of an efficient 
and unbiased estimator (i.e., an estimation method that achieves the CRB on parameter 
precision). This has allowed us to calculate lower limits on the impact of point source 
subtraction, and using a first principles approach, but does not address the issue of the 
existence of such a method (a method may not exist that can achieve the bound). In 
addition to designing an efficient estimator, the CRB applies only to unbiased estimators, 
and introducing bias into the estimation will yield different structures in the final power 
spectrum than shown here. As a simple estimate of the impact of an unbiased, but inefficient 
estimator, we (arbitrarily) assume we have estimated the source parameters to within three 
times the CRB, and calculate the residual signal uncertainty in the angular power spectrum. 
This toy model yields an increase in uncertainty of a factor of ~10, indicating a rough square 
scaling of estimator inefficiency with noise uncertainty. 

We have derived a ~|w| 2 -dependence of visibility-based variance on baseline length, 
demonstrating a ramping of noise power at large angular modes (small angular scales). In 
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this work, we have used the same antennas to estimate the source parameters and produce the 
power spectra (the only difference being the bandwidths used). The long baselines contribute 



relatively more than the short baselines to the estimate of source position (equations 13-14). 
However, these baselines form the largest I modes, where residual signal has the greatest 
impact. Inclusion of these baselines for parameter estimation is crucial for precise estimates; 
their use in EoR measurements depends on the experimental design, and the desire for 
measuring the 21 cm EoR signal at those scales. Use of these baselines does not affect EoR 
estimation at other scales (the correlation between modes is small). Also, as demonstrated 
Trott et al. ( 2011[ ), inclusion of short baselines for position estimation is important (but 



m 



less so than long baselines), due primarily to the increased instrumental sensitivity when 
using all available antennas. 

Our results also demonstrate the expected behaviour of minimally-redundant versus 
maximally-redundant arrays. Uniform sampling yields even thermal noise power across a 
wide range of angular modes, while maximally-redundant arrays sample some scales well, 
but yield patchy performance and restricted range of sampled modes. Choice between these 
models depends on the experiment being performed, but a maximally-redundant array re- 
quires prior information about the most interesting scales for estimating the 21 cm signal, 
while the uniform arrays are a better "starter" design, where the scientist has no prior 
information about signal structure and is aiming for a first detection/estimation. 



Appendix 

A. Scaling of expected power spectrum and variance with number of 

visibilities 

Computation of the angular power spectrum involves both a coherent and incoherent 
addition of visibilities. The relative size of the coherent to incoherent calculation depends 
on the number of visibilities contributing to each (uu)-cell, and the number of (uv) cells 
contributing to each /-mode. In general, the relative size of these additions is different 
for each /-mode. As such, the scaling of the expected value of the point source residual 
power spectrum, and its variance, with the number of contributing visibilities is non-trivial, 
and the design of the interferometer, and the binning of the power spectrum, affect the 
result. In this Appendix, we will briefly outline how the relative size of the coherent to 
incoherent combination scales with the number of contributing visibilities, beginning with 
the two extreme examples (all coherent and all incoherent), and then for an intermediate 
example. This is motivated by the discussion of minimum-redundancy versus maximum- 



redundancy configurations of Parsons et al. (2011), and the initial use of the PAPER array 
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in both types of configurations. Typically, visibilities can be considered coherent over scales 
that approximate the instrument field-of-view. Arrays that sample the same angular modes 
coherently (maximally-redundant) can reduce the contribution of thermal noise on these 
scales, at the expense of being able to sample a wider range of angular scales (for a fixed 
number of antennas). 

The angular power spectrum is estimated using a maximum-likelihood estimator, de- 
signed to weight the contributions from each cell according to the noise they contain. The 
estimate for mode / is given by; 



|2 



C = ^-. (AD 

/ uv 

(uv)ei 

where N uv is the number of visibilities contributing to a given (uv) cell, and the sum is over 
the (uv)-ceWs contributing to that /-mode. The visibilities, V uv , are coherently combined 
(averaged) within each (uv) cell, while the final estimate associating (m>)-cells in /-modes 
is incoherent (visibilities are combined after squaring). For simplicity, we will compute the 
statistics for complex white Gaussian noise on each visibility, distributed as V ~ £/V(0, cr 2 ) 
(i.e., the thermal noise, where l CM ~ (0, cr 2 )' denotes that the PDF is distributed according 
to a zero mean complex Gaussian distribution with variance, a 2 ). 

For an entirely coherent estimate (jVj = N uv ), where a single (uv) cell contributes to a 
single /-mode, the summations are omitted; 

Q = \V UV \ 2 , (A2) 

yielding, 



C>~« r^)- (A3) 



V 2a^ 
Nl 'N 2 

Note that in this case, increasing the number of visibilities reduces the expected noise power 
and variance, but the ratio of expected noise power to standard deviation of power is un- 
changed. Strictly speaking, this is a chi-square distribution with one degree of freedom, 
yielding a highly skewed probability distribution function, /, 

m) ^ » (A4) 



where x = C\o 2 /N t . 
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At the opposite extreme, for an entirely incoherent estimate (N uv = 1), where there is 
a single visibility contributing to each (uv) cell, the power spectrum estimate is given by: 

I>l 2 



yielding, 



C = ^— , (A5) 



Q ~ M ■ (A6) 



In the intermediate case, where we would expect to be for real experiments, there will 
be large numbers of visibilities combined coherently, and a mapping of lOs-lOOs (uv) cells to 
each I mode. We take an example with ten (uv) cells contributing to a given /-mode, with 
an even number of visibilities in each (w)-cell, Ni = 10N UV . The power spectrum estimate 
is given by: 



10 

|2 



^N UV \V W 

d = ^ , (A7) 

V N 

/ uv 



C~Af(#,^). (AS) 



yielding, 

a ~ .am 

This expression forms the natural transition from completely coherent to completely inco- 
herent, and we have again taken the Gaussian approximation to the chi-square distribution. 



B. Pedagogical motivation for the 'wedge' feature 



The 'wedge' of residual power expected from point source subtraction has been discussed 



by Datta et al. (2010), Vedantham et al. (2012), and Morales et al. (2012), and represents the 



power footprint of the instrument point spread function. In particular, the review of Morales 



et al. (2012) presents arguments for the origin of this feature due to point source amplitude 



and position errors, gridding artifacts and calibration errors. They motivate the position 
and structure of the wedge, but without a first-principles approach to the magnitude of the 
errors. In this Appendix, we extend the previous discussion by (1) motivating the shape of the 
wedge from the first-principles approach followed here, providing an alternative perspective 
on its shape and location, and (2) complementing previous descriptions by describing the 
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integrated effect of multiple imprecisely-subtracted point sources. Note that we are not 
suggesting previous analyses are incorrect; we are simply approaching the problem from a 
different path. 

We have already demonstrated that with an unbiased and efficient estimate of the 
parameters of each point source (i.e., obtaining the CRB on parameter precision), the error 
introduced into a given visibility scales as the square of the baseline length (see Equations 



23 in Section 3.3). This analysis explains the ramping of the residual signal power in the 
k± dimension. In this Appendix we focus on the shape in the k\\ direction; specifically, the 
sharp drop-off of signal from small to large fcn at large angular scales (small k±), and the 
relatively flat signal power at small angular scales (large k±). We will demonstrate that this 
structure is a combination of correlated signal across frequency channels spreading across 
multiple uv -cells ("mode- mixing"), and the integrated effect of subtracting multiple sources 
distributed randomly across the sky. 

For a single source contaminating a (u, v, v) cell, the magnitude of the covariance be- 
tween frequency channels (for the same (u,v)) is equal to unity (i.e., they are completely 
correlated). There is, however, a phase shift in the complex covariance matrix, yielding the 



following structure (equation 20): 

/ 1 

exp 



4 



exp 



—%<p 
-2k 



expzc; 
1 

exp —i 



exp 2i<j 
exp iff) 
1 



\exp-Ni(j) exp-(iV - l)i<j) exp -(N - 2)i<j) 



exp Ni(f> \ 
exp(iV - l)i(j> 
exp (N - 2)i0 



1 



(Bl) 



/ 



where 



= 2vrAz/(AxZ + Aym)/c (B2) 

is the phase term over frequency range Au (a single channel) for a source located at sky 
position (l,m) and baseline lengths (Ax, Ay), and aj> is the source-independent variance in 



a visibility, given by equation [23] This equation can be derived using the off-diagonal entries 
of equation [20] 



When there are multiple sources contributing residual signal to the coherently-averaged 
visibility, these phase terms add, yielding a sum over independent phases for the matrix 
entry describing the covariance between frequency channels, u a ,vp: 

C uv (u a , up) = J^exp (a - fiifa. (B3) 

3=1 

For phases distributed evenly over 2tc, the expected value of this sum is zero, and the 
frequency channels would be completely uncorrelated. However, the visibility phase term, (p, 
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is limited in magnitude by the extent of the array (maximum (Ax, Ay)) and the field-of-view 
(maximum (/, to)). If we assume, for simplification, that the phases are uniformly-distributed 
over some range, [— 9, 9} (corresponding to randomly-located sources and a uniform density 
of antennas), the expected value of the summed phases is: 



where 9 = max , and the variances (matrix diagonal entries) have also summed to iV ca i. 
Note that the expected value of the off-diagonal terms sum to a real quantity. Hence, 
the covariance matrix, summed over multiple sources, is real, with a correlation between 
frequency channels that decays as a sine function with an argument that is linear in the 
frequency spacing between the channels under consideration. The maximum value of is 
determined by the instrument field-of-view and the baseline length under consideration. At 
small values of \u\ 2 (i.e., at small k± values) the maximum value is smaller than at larger 
baseline lengths (large k± values). Hence, the correlation between frequency channels decays 
more rapidly at small angular scales (large k±). For the MWA, with a maximum baseline 
of 3000 m (minimum 5 m), a field-of-view of 30°, and frequency channels of 80 kHz (used in 
this work), the maximum values of for the shortest and longest baseline are ~0.002 and 
~1.3, respectively. Therefore, on the longest baselines, the sine function reaches its first null 
after a few channels, while on the shortest baselines, it is >1000 (which exceeds the number 
of channels considered for a 8 MHz bandwidth). 

In addition to the integrated effect of imprecisely subtracting point sources, the noise 
structure of the power spectrum is affected by mode-mixing: the evolution of the instrument 
PSF with frequency, which manifests itself in distributing the signal from a given baseline 
over several (uv) cells. In the context of covariances between frequency channels in the same 
(uv ) cell, this straying of visibilities leads to a cutoff in the correlation length, dependent 
upon the width of the coherent (uv ) cell relative to gradient in (uv) evolution with frequency. 
Mode-mixing is negligible on short baselines, and can be considerable on long baselines. The 
size of a coherent (uv ) cell is given approximately by the inverse of the instrument field-of- 
view, yielding a Aw cc n ~ 2 cell-size for the MWA with a 30° field. The spread of u values 
for a given baseline is given by Au = ArrBW/c, which has the range [0.1,80] for the MWA. 
Hence, at the shortest baselines, all visibilities are contained within a single cell, while they 
begin to shift cells for Ax > Aw cell c/BW « 70 m (with the longest baselines traversing >40 
cells). 

Armed with an understanding of the structure of the frequency-based covariance matrix, 
we can now explore how the covariances propagate through the v — r] Fourier transform and 
power (mod-squared) operations, and lead to the wedge feature observed in the results. For 



E 




(B4) 
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an uncorrelated dataset, the variance in the power formed via DFT and squaring is the 
same for all Fourier modes (77), even for different variances in the individual frequency-based 
dataset, and is given by the sum of the squares of the variances of each component (equation 



24]: 

^ = 2tr(ACUC) = (B5) 

i=i 

where N c h = N v is the number of frequency channels and Fourier bins. This behaviour 
corresponds to the long baselines, where the combination of mode-mixing and small frequency 
correlation length yield relatively uncorrelated matrices. Visibilities on the longest baselines 
traverse cells after ~ 2 — 3 frequency channels, yielding a substantially truncated channel 
correlation, and effectively uncorrelated frequency channels. At the other extreme, the larger 
angular modes (small k±, uv) suffer from little mode-mixing and have longer correlation 
lengths, yielding qualitatively different behaviour in the DFT+squaring operations. 

We now motivate the shape of the wedge as a function of k±. We have demonstrated that 
the frequency channels are correlated with a sine function dependence, yielding a visibility 
dataset that corresponds to Gaussian noise convolved with a sine function. The Fourier 
transform of this convolution yields a flat line-of-sight power, multiplied by the Fourier 
transform of a sine function (a rectangle function). I.e., it performs a low-pass filtering 
operation. Hence, the expected noise power follows a rectangle function (squared), with a 
width determined by the correlation length in the frequency domain: 

F[M(0, a 2 ) * sinc(0 max n)] (B6) 
= F[AT(O,a 2 )].F[sinc(0 max n)] (B7) 

oc F[sinc(0 max n)] (B8) 

f nk \ . . 

oc rect — , (B9) 



max 



where 

rect(x) = {' ----- ( B10 ) 




is the rectangle function, and n,k = [0,1,.., AT — 1] index the DFT. This function yields 
a constant expected power when k = N, and a single bin rectangle function when k — 1, 



producing the following limits on baseline lengths (using equation B2): 



Ax < 70m if jfe = 1 (Bll) 
Ax ~ 7000m if jfe = N. (B12) 
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The latter baseline length exceeds the maximum for the MWA. One would therefore expect 
that the noise power would not be flat across the line-of-sight range, but instead fall off 
at large k\\. This analysis neglects the mode-mixing, however, and the sine function (and 
therefore the correlation length) is severely truncated for the longest baselines. The sine 
function is always truncated at its first null (argument=7r), because this corresponds to the 
scale over which the visibilities for a given baseline traverse the coherent (uv) cell. This can 
be seen by considering the number of frequency channels to reach the first sine function null, 
and the number of frequency channels over which a baseline is contained within a single cell. 
In one dimension: 



-^Va lysine 

Na 

1 ' AzAmm 



Au r 



2Aur ma , x Ax 
1 c 



2r n 



Ai> Ax ' 



(B13) 
(B14) 



where 'mm' refers to mode- mixing. These two expressions are seen to be equivalent. The 
truncation windows the frequency-domain sine correlation with a rectangle function, tapering 
the power-domain rectangle function with a sinc-like function: 



oc 



T [sinc(n0)] 
nk 

1; ° C1 ' ^0 

7T 2 k 



rect 



•k sine 



ncj) 

IT 
7l 2 k 
N^4> 



~ sine 



(B15) 
(B16) 
(B17) 



where k denotes the line-of-sight Fourier bin. Inserting the expression for = 4> max and 
transforming co-ordinates to k± and k\\ yields the following approximate expression for the 
wedge structure: 

^^||) °c sine 2 f — (B18) 



k± BW 



where, 



ttc(1 + z) 2 D M (z) 



(B19) 

are array-independent terms defining the cosmology and observation redshift (with distance 
units), BW is the experiment bandwidth at observation frequency, v. Here W represents 
the uncertainty in noise power (the expected power follows the same functional form) due 
to imprecise subtraction of multiple point sources, spread randomly across the instrument 
field-of-view. 



We have therefore arrived at the final form of the frequency covariance matrix — a 
real-valued truncated sine function, with a truncation length at the first null — yielding 
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a sinc-squared-like function with characteristic length given by the argument of equation 
B18| Figure [9] plots examples of this function, as a function of line-of-sight Fourier bin, for 
different values of k±. The analytic result mimics the qualitative structure observed in the 

Uncertainty in noise power 




1 00 



Log k bin 



Fig. 9. — Uncertainty in noise power for different values of k± for the MWA (normalized to 
first non-DC term in the power). 

power spectra shown in this work, and in previous published results (the result is approximate 
because it relies on simplifying assumptions, namely that the correlation length is truncated 
exactly at the first null of the sine function, and that the antenna layout is such to produce 
a uniform sampling of correlation phases in the frequency covariance matrix). 
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